{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "26c36742-9918-4270-8f11-b6fa23eb36f5",
   "metadata": {},
   "source": [
    "<img src=\"../common/rfsoc_book_banner.jpg\" alt=\"University of Strathclyde\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d993d47-165c-4ac2-a488-6adcf779337b",
   "metadata": {},
   "source": [
    "# Notebook Set C\n",
    "\n",
    "---\n",
    "\n",
    "## 03 - Fourier's Theorem\n",
    "In this notebook, we will investigate Fourier's Theorem and the Complex Fourier Series. We will begin by investigating a sum of sine and cosine waves to approximate a square wave. Then, we will define a function for generating the Fourier Series for any real arbitrary waveform. This function will be used to determine the Fourier Series of a sawtooth wave. Euler's formula will also be introduced and a function for a Complex Fourier Series will be defined. Lastly, we will explore a simple example of obtaining the Complex Fourier Series coefficients for an arbitrary waveform. \n",
    "\n",
    "## Table of Contents\n",
    "* [1. Introduction](#introduction)\n",
    "* [2. A Sum of Sine and Cosine Waves](#a_sum_of_sine_and_cosine_waves)\n",
    "* [3. Defining a Fourier Series Function](#defining_a_fourier_series_function)\n",
    "* [4. Euler's Formula](#eulers_formula)\n",
    "* [5. The Complex Fourier Series](#complex_fourier_series)\n",
    "* [6. Conclusion](#conclusion)\n",
    "\n",
    "    \n",
    "## References\n",
    "* [1] - [J. Fourier, “Théorie analytique de la chaleur,” (in French) Paris: Firmin Didot Père et Fils, 1822.](https://archive.org/details/bub_gb_TDQJAAAAIAAJ)\n",
    "* [2] - [E. Hewitt and R. E. Hewitt, “The Gibbs-Wilbraham Phenomenon: An Episode in Fourier Analysis,” Archive for History of Exact Sciences, Vol. 21, 1979, pp. 129-160.](http://doi.org//10.1007/BF00330404)\n",
    "\n",
    "## Revision History\n",
    "* **v1.0** | 21/12/2022 | _First Revision_\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8b2669c-2f05-4bff-9dd1-23dd75804d26",
   "metadata": {},
   "source": [
    "## 1. Introduction <a class=\"anchor\" id=\"introduction\"></a>\n",
    "Joseph Fourier presented the process for decomposing a wave into a set of harmonically related periodic waveforms in his book titled La Théorie Analitique de la Chaleur [1]. This decomposition is known as a Fourier Series. Any arbitrary periodic waveform of fundamental period $T$ can be decomposed into a sum of sine and cosine waves of varying amplitudes and harmonically related frequencies. The following code cells will demonstrate how to obtain the Fourier Series of any arbitrary waveform. By the end of this notebook, you will have a greater understanding of practical Fourier Series analysis.\n",
    "\n",
    "## 2. A Sum of Sine and Cosine Waves <a class=\"anchor\" id=\"a_sum_of_sine_and_cosine_waves\"></a>\n",
    "The Fourier Series of a periodic waveform is defined as\n",
    "\n",
    "$$\n",
    "g(x) = a_0 + \\sum_{n=1}^{\\infty} a_n \\cos{(nx)} + \\sum_{n=1}^{\\infty} b_n \\sin{(nx)}.\n",
    "$$\n",
    "\n",
    "Let us begin our investigation by obtaining the Fourier Series of a square wave. Start by importing the required libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68bddf66-3a4b-4db2-8366-87e5970bf9cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import signal, integrate\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1da0a95-524c-4ee3-8892-2ffe1d122627",
   "metadata": {},
   "source": [
    "Now, define a function for a square wave and plot it for inspection. The square wave will be characterised by the following:\n",
    "* sample frequency $f_s$ of 1kHz,\n",
    "* a fundamental period $T$ of 2$\\pi$,\n",
    "* an amplitude of 3,\n",
    "* a DC bias of +3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13657eda-3b21-45eb-9ef8-b3449269cf27",
   "metadata": {},
   "outputs": [],
   "source": [
    "def squarewave(x, amplitude=3, dcbias=3, phase=np.pi/4):\n",
    "    return signal.square(x+phase)*amplitude+dcbias\n",
    "\n",
    "fs = 1e3\n",
    "x = np.arange(-2*np.pi, 2*np.pi, 1/fs)\n",
    "plt.plot(x, squarewave(x), label='Square Wave', linestyle='dashed')\n",
    "plt.ylim([-1, 7])\n",
    "plt.xlim([-2*np.pi, 2*np.pi])\n",
    "plt.xlabel('Phase (Radians)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbb9ade8-729d-4076-9aaf-2a6ecafc51f6",
   "metadata": {},
   "source": [
    "The first Fourier Series coefficient, denoted as $a_0$, is the DC bias of a periodic waveform and is calculated using\n",
    "\n",
    "$$\n",
    "a_0 = \\frac{1}{2\\pi}\\int\\limits_{-\\pi}^\\pi g(x) \\,dx.\n",
    "$$\n",
    "\n",
    "We can leverage the scipy `integrate.quad` function to compute $a_0$ across the periodic interval $[-\\pi, \\pi]$ and obtain the DC bias of the square wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e0fe3d5-e3ec-4a6d-93db-1eef792f5b9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "a0 = integrate.quad(squarewave, -np.pi, np.pi)[0]/(2*np.pi)\n",
    "print('The DC bias is ' + str(a0) + '.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a628ee08-14ae-4bc3-a17e-5fea43c1d4fb",
   "metadata": {},
   "source": [
    "The cosine and sine Fourier Series coefficients are denoted as $a_n$ and $b_n$, respectively, where $n=1,2,3, ...$ The cosine coefficients are defined as\n",
    "\n",
    "$$\n",
    "a_n = \\frac{1}{\\pi}\\int\\limits_{-\\pi}^\\pi g(x)\\cos{(nx)} \\,dx.\n",
    "$$\n",
    "\n",
    "The sine coefficients are computed using\n",
    "\n",
    "$$\n",
    "b_n = \\frac{1}{\\pi}\\int\\limits_{-\\pi}^\\pi g(x)\\sin{(nx)} \\,dx.\n",
    "$$\n",
    "\n",
    "We can define functions for each of the terms inside the integrals of the above equations. Then, integrate the functions over the interval $[-\\pi, \\pi]$ to obtain the first Fourier Series coefficients for the cosine and sine waves, which are denoted as $a_1$ and $b_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4a90569-cc95-4355-a53c-bb936a977f34",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gcos(x, g, n):\n",
    "    return g(x)*np.cos(n*x)\n",
    "\n",
    "def gsin(x, g, n):\n",
    "    return g(x)*np.sin(n*x)\n",
    "\n",
    "g = squarewave\n",
    "n = 1\n",
    "a1 = integrate.quad(gcos, -np.pi, np.pi, args=(g, n))[0]/(np.pi)\n",
    "b1 = integrate.quad(gsin, -np.pi, np.pi, args=(g, n))[0]/(np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb4a48c9-0bc8-4e91-a999-d4a4202c3c86",
   "metadata": {},
   "source": [
    "The DC bias, cosine wave, and sine wave can be summed together as below. The result is then plotted alongside the square wave for inspection. We can clearly see a waveform that approximately matches the amplitude and periodicity of the square wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afb45c46-0d1f-45cb-bbb4-fff004094664",
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = a0 + a1*np.cos(n*x) + b1*np.sin(n*x)\n",
    "plt.plot(x, squarewave(x), label='Square Wave', linestyle='dashed')\n",
    "plt.plot(x, s1, label='1 Term')\n",
    "plt.ylim([-1, 7])\n",
    "plt.xlim([-2*np.pi, 2*np.pi])\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('Phase (Radians)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84772dd2-0e53-4342-9aab-d7f271ae291d",
   "metadata": {},
   "source": [
    "To obtain a better approximation, we can increase the number of Fourier Series terms that are summed together to 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a90fbf52-4c6c-4672-9861-f8314e608958",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 3\n",
    "a3 = integrate.quad(gcos, -np.pi, np.pi, args=(g, n))[0]/np.pi\n",
    "b3 = integrate.quad(gsin, -np.pi, np.pi, args=(g, n))[0]/np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e333742-36b7-4028-809c-ab07505a075c",
   "metadata": {},
   "source": [
    "Plotting the summation of 3 terms improves our square wave approximation, as given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57536e3d-1368-4b9e-842a-34bf5e12cfe3",
   "metadata": {},
   "outputs": [],
   "source": [
    "s3 = s1 + a3*np.cos(n*x) + b3*np.sin(n*x)\n",
    "plt.plot(x, squarewave(x), label='Square Wave', linestyle='dashed')\n",
    "plt.plot(x, s3, label='3 Terms')\n",
    "plt.ylim([-1, 7])\n",
    "plt.xlim([-2*np.pi, 2*np.pi])\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('Phase (Radians)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69ae62f4-95b6-4c1a-900e-dfb608c150ec",
   "metadata": {},
   "source": [
    "A function can be defined to increase the number of Fourier Series terms that we sum across. We create this function in the following section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c01a29d-ec3d-4850-9fc1-e60571ab5645",
   "metadata": {},
   "source": [
    "## 3. Defining a Fourier Series Function <a class=\"anchor\" id=\"defining_a_fourier_series_function\"></a>\n",
    "We can define a function to compute the Fourier Series coefficients for any number of terms and across any periodic interval. This function will compute the Fourier Series for any arbitrary real valued waveform and is given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "223b8919-9121-4d0f-b157-d52ac2ee1926",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fourier_coeffs(x, g, n, a=-np.pi, b=np.pi):\n",
    "    an, bn = np.zeros(n), np.zeros(n)\n",
    "    a0 = integrate.quad(g, a, b)[0]/(b-a)\n",
    "    for i in range(1, n + 1):\n",
    "        an[i-1] = integrate.quad(gcos, a, b, args=(g, i))[0]*2/(b-a)\n",
    "        bn[i-1] = integrate.quad(gsin, a, b, args=(g, i))[0]*2/(b-a)\n",
    "    return a0, an, bn"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "736e1415-3874-4890-8bcc-0b472f0873c6",
   "metadata": {},
   "source": [
    "Once we have calculated the Fourier Series coefficients, we can sum them together to produce the Fourier Series approximation of the waveform. We can define another function to perform the summation as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5ffd81d-24dd-4892-bff7-3994d08e4aa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fourier_sum(x, a0, an, bn, n=None):\n",
    "    if n == None or n > len(an):\n",
    "        n = len(an)\n",
    "    s = np.ones(len(x))*a0\n",
    "    for i in range(1, n + 1):\n",
    "        s += an[i-1]*np.cos(i*x) + bn[i-1]*np.sin(i*x)\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15922fc1-1390-4f9c-81e0-1096adcf7c8c",
   "metadata": {},
   "source": [
    "The functions can now be used to obtain the Fourier Series of the square wave for any number of terms. Let us begin by computing the first 26 Fourier Series coefficients using the `fourier_coeffs` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6072b23-663e-48b7-92a9-d49f17dc1915",
   "metadata": {},
   "outputs": [],
   "source": [
    "a0, an, bn = fourier_coeffs(x, squarewave, n=26)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa68b9c1-c58f-4b73-a123-40f8431854d8",
   "metadata": {},
   "source": [
    "Now, we can compute the sum of the first 13 terms and plot the result as given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "606b2627-5d20-4b5b-8060-4c0dafd24717",
   "metadata": {},
   "outputs": [],
   "source": [
    "s = fourier_sum(x, a0, an, bn, n=13)\n",
    "plt.plot(x, squarewave(x), label='Square Wave', linestyle='dashed')\n",
    "plt.plot(x, s, label='13 Terms')\n",
    "plt.ylim([-1, 7])\n",
    "plt.xlim([-2*np.pi, 2*np.pi])\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('Phase (Radians)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17f075a3-9214-4b61-8bd4-d99e1f5cd5a5",
   "metadata": {},
   "source": [
    "To improve the approximation, we can sum all 26 terms of the Fourier Series expansion. Notice that the resulting waveform is closely resembling the square wave. The ringing near the square wave transitions is known as Gibbs Phenomenon, which can be investigated further in [2]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "409de604-c678-4c60-a2e1-494efdbc9bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "s = fourier_sum(x, a0, an, bn, n=26)\n",
    "plt.plot(x, squarewave(x), label='Square Wave', linestyle='dashed')\n",
    "plt.plot(x, s, label='26 Terms')\n",
    "plt.ylim([-1, 7])\n",
    "plt.xlim([-2*np.pi, 2*np.pi])\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('Phase (Radians)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d648cf6-67f7-418f-90f4-3e1272cc9b7e",
   "metadata": {},
   "source": [
    "These functions can also be used to obtain the Fourier Series of other waveforms, such as a sawtooth wave. Firstly, let us define a function for a sawtooth waveform, as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7f04b2f-7267-4bc5-9d60-765c04b4996a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sawtoothwave(x, amplitude=0.5, dcbias=0.5, phase=np.pi):\n",
    "    return signal.sawtooth(x+phase)*amplitude+dcbias"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8db4f19-af82-4001-88eb-ea38365c978b",
   "metadata": {},
   "source": [
    "Now we can create a sawtooth waveform that features the following properties:\n",
    "* sample frequency $f_s$ of 1kHz,\n",
    "* a fundamental period $T$ of 1 across the interval [0, 1] in seconds,\n",
    "* an amplitude of 0.5,\n",
    "* a DC bias of +0.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4068d19b-6eb5-498f-ac07-108d1c9309f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 1e3\n",
    "T = 1/2\n",
    "ts = np.arange(0, 1, 1/fs)\n",
    "x = 2*np.pi*ts/T\n",
    "plt.plot(ts, sawtoothwave(x), label='Sawtooth Wave', linestyle='dashed')\n",
    "plt.ylim([-0.2, 1.2])\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8b60d54-3c2b-4c38-87eb-69d41bdda432",
   "metadata": {},
   "source": [
    "We can use the `fourier_coeffs` function to compute the first 5 terms of the fourier series expansion. Also, the `fourier_sum` function can be used to obtain the sum after computing each term. We can plot the summations on the same graph as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "decfabba-0402-4c57-a7ac-0adc95f5f822",
   "metadata": {},
   "outputs": [],
   "source": [
    "a0, an, bn = fourier_coeffs(x, sawtoothwave, n=5, a=x[0], b=x[-1])\n",
    "s1 = fourier_sum(x, a0, an, bn, n=1)\n",
    "s2 = fourier_sum(x, a0, an, bn, n=2)\n",
    "s3 = fourier_sum(x, a0, an, bn, n=3)\n",
    "plt.plot(ts, sawtoothwave(x), label='Sawtooth Wave', linestyle='dashed')\n",
    "plt.plot(ts, s1, label='1 Term')\n",
    "plt.plot(ts, s2, label='2 Terms')\n",
    "plt.plot(ts, s3, label='3 Terms')\n",
    "plt.ylim([-0.2, 1.2])\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3d3d46c-13cf-46f3-ba49-393ddf26441c",
   "metadata": {},
   "source": [
    "In this section, we have obtained the Fourier Series of a square and sawtooth waveforms. We have seen that a sum of harmonically related sine and cosine waves can be used to approximate any arbitrary real waveform.\n",
    "\n",
    "We will soon introduce the Complex Fourier Series to compute the Fourier Series coefficients of a complex waveform. Before we investigate the Complex Fourier Series, we will first explore Euler's formula in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16643019-8b39-4f93-b209-cfb8452a8d1b",
   "metadata": {},
   "source": [
    "## 4. Euler's Formula <a class=\"anchor\" id=\"eulers_formula\"></a>\n",
    "Euler's formula describes the relationship between a complex exponential function and associated trigonometric functions. Euler's formula is commonly used to simplify the mathematical notation of complex numbers and is defined as,\n",
    "\n",
    "$$\n",
    "e^{j\\theta} = \\cos{(\\theta)} + j\\sin{(\\theta)}.\n",
    "$$\n",
    "\n",
    "The relationship uses the base of natural logarithms $e$, which is also known as Euler's constant. The numpy library defines Euler's constant as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a050c09c-30b6-4da9-8ccb-2d51003578ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.e"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c792d19-0e87-44d5-8c09-6bc290c701ef",
   "metadata": {},
   "source": [
    "We can define a complex number using the trigonometric relationship below, where the angle is $\\pi/8$ and the amplitude is 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35955ed8-f047-46c6-8d1d-fcc3042a1de4",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.cos(np.pi/8) + 1j*np.sin(np.pi/8)\n",
    "np.round(c, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a174f7ef-dc02-4572-b0c0-8bb50e8f6b62",
   "metadata": {},
   "source": [
    "In a similar way, we can also use Euler's constant to define the same complex number. Running the next code cell demonstrates that Euler's relationship holds true for these values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45c07931-b352-4464-b6d3-06f84f8e4e0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.e**(1j*np.pi/8)\n",
    "np.round(c, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e62a8fc-87fc-4944-a6ee-93a641f7e896",
   "metadata": {},
   "source": [
    "Euler's formula can also be used to compute the complex conjugate, as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a70e4225-e64e-49e2-bf47-21da784643e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "cbar = np.e**(-1j*np.pi/8)\n",
    "np.round(cbar, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6f43255-c1c6-4baa-849a-b74b21a51f6b",
   "metadata": {},
   "source": [
    "Euler's formula is very useful as it can reduce the difficulty of performing mathematical operations on complex numbers. Multiplying complex exponentials is simple, as we only need to sum their powers. For instance,\n",
    "\n",
    "$$\n",
    "e^{(A + B)} = e^{A} e^{B}\n",
    "$$\n",
    "\n",
    "An example is given in the following code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c4f1350-081f-4f5f-86e0-a13356970949",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.e**(2) * np.e**(3) == np.e**(2 + 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a0662b0-1af4-4fa0-a679-b54ce24bfaa8",
   "metadata": {},
   "source": [
    "Now that we have an understanding of Euler's formula, we can progress to the Complex Fourier Series."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff9687d9-2ba6-43e0-8f11-77aa4a643a55",
   "metadata": {},
   "source": [
    "## 5. The Complex Fourier Series <a class=\"anchor\" id=\"the_complex_fourier_series\"></a>\n",
    "Periodic complex waveforms can also be decomposed into a sum of harmonically related sine and cosine waves. The generalised Complex Fourier Series for an arbitrary complex waveform can be defined as,\n",
    "\n",
    "$$\n",
    "g(t) = \\sum_{n=-\\infty}^{\\infty} c_n e^{j2\\pi nt/T},\n",
    "$$\n",
    "\n",
    "where $n = 0, \\pm 1, \\pm 2, ...$ Let us begin by defining a function for a complex waveform. For example, we can set the real part to a sawtooth wave and the imaginary part to a square wave, as given in the function below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9209e6c3-d345-4e52-9e89-912ee23ac87c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def complex_wave(x, amplitude=1+0.5j, dcbias=0+0j, phase=np.pi/4+0j):\n",
    "    re = signal.sawtooth(x+np.real(phase)+np.pi)*np.real(amplitude)+np.real(dcbias)\n",
    "    im = signal.square(x+np.imag(phase))*np.imag(amplitude)+np.imag(dcbias)\n",
    "    return re + 1j*im"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1f47855-ac4b-4645-9d9a-736e9f75827b",
   "metadata": {},
   "source": [
    "We can define a complex waveform using the above function and plot it for inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a3dd98a-bdaa-478e-90be-5ee09540b547",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 1e3\n",
    "ts = np.arange(0, 0.5, 1/fs)\n",
    "T = 1/4\n",
    "x = 2*np.pi*ts/T\n",
    "plt.plot(ts, np.real(complex_wave(x)), label='Real Wave', linestyle='dashed')\n",
    "plt.plot(ts, np.imag(complex_wave(x)), label='Imag Wave')\n",
    "plt.ylim([-1.2, 1.2])\n",
    "plt.xlim([0.0, 0.5])\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "925146b7-de71-401c-892b-0688de0d3e7c",
   "metadata": {},
   "source": [
    "The Complex Fourier Series coefficients are denoted as $c_n$, which can be computed using the integral,\n",
    "\n",
    "$$\n",
    "c_n = \\frac{1}{T} \\int_\\limits{-\\infty}^\\infty g(t)e^{-j2\\pi f_0 nt} \\, dt.\n",
    "$$\n",
    "\n",
    "The term inside the integral above can be calculated using the function below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "613fa8f3-1dca-4479-bba2-8e860c6cf97f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gexp(x, g, n):\n",
    "    return g(x)*np.e**(-1j*n*x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dc2a885-88c0-498c-90ab-96d8b1250810",
   "metadata": {},
   "source": [
    "As before for the Fourier Series, we can create functions to compute the Complex Fourier Series coefficients and summation. Also, we have introduced functions that extract the real and imaginary components from the `gexp` function above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd17bac7-deeb-4aef-84b1-ed722c76dda7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def real_gexp(x, g, n):\n",
    "    return np.real(gexp(x, g, n))\n",
    "\n",
    "def imag_gexp(x, g, n):\n",
    "    return np.imag(gexp(x, g, n))\n",
    "\n",
    "def complex_fourier_coeffs(x, g, n, a=-np.pi, b=np.pi):\n",
    "    an, bn = np.zeros(2*n+1), np.zeros(2*n+1)\n",
    "    for i in range(-n, n+1):\n",
    "        an[i+n]=integrate.quad(real_gexp, a, b, args=(g, i))[0]/(b-a)\n",
    "        bn[i+n]=integrate.quad(imag_gexp, a, b, args=(g, i))[0]/(b-a)\n",
    "    cn = an + 1j*bn\n",
    "    return cn\n",
    "\n",
    "def complex_fourier_sum(x, cn):\n",
    "    n = int(np.ceil(len(cn)/2))-1\n",
    "    s = np.zeros(len(x), dtype=complex)\n",
    "    for i in range(-n, n+1):\n",
    "        s += (cn[i+n]) * np.e**(1j*x*i)\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e26c75ff-d0c3-4680-b139-1d057ac2f09f",
   "metadata": {},
   "source": [
    "Now, we can compute the Complex Fourier Series coefficients for the complex waveform. The code cell below computes the first 9 terms and plots the final summation for inspection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fceb76f3-0366-4a88-aa55-e0dede557ae8",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 9\n",
    "cn = complex_fourier_coeffs(x, complex_wave, n)\n",
    "s = complex_fourier_sum(x, cn)\n",
    "plt.plot(ts, np.real(s), label='Real Wave (9 Terms)', linestyle='dashed')\n",
    "plt.plot(ts, np.imag(s), label='Imag Wave (9 Terms)')\n",
    "plt.ylim([-1.2, 1.2])\n",
    "plt.xlim([0.0, 0.5])\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "655dabdd-32ee-43ba-945e-9de3501369f5",
   "metadata": {},
   "source": [
    "You can clearly see the sawtooth wave and square wave approximations. The Complex Fourier Series reduces the difficulty of decomposing complex waveforms into a sum of harmonically related sine and cosine waves. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36901568-644a-4073-9ead-3f7ae7224a3e",
   "metadata": {},
   "source": [
    "## 6. Conclusion <a class=\"anchor\" id=\"conclusion\"></a>\n",
    "In this notebook, we investigated the Fourier Series and the Complex Fourier Series for the harmonic analysis of periodic waveforms. In the following notebook, we will explore the Discrete Fourier Transform (DFT) and its operating principles."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "250aa0fa-a7de-4fb3-a080-7f5827cf36fa",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "[⬅️ Previous Notebook](02_rfsoc_spectrum_analyser.ipynb) || [Next Notebook 🚀](04_investigating_the_dft.ipynb)\n",
    "\n",
    "Copyright © 2023 Strathclyde Academic Media\n",
    "\n",
    "---\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
